
/********************************** LICENSE ************************************
 *
 * Copyright 2011 Alessandro Carrega
 * email: alessandro.carrega@unige.it
 *
 * This file is part of model_param_simulator.
 *
 * model_param_simulator is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * model_param_simulator is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with model_param_simulator.  If not, see <http://www.gnu.org/licenses/>.
 *
 ********************************* LICENSE ************************************/

#include <ginac/ginac.h>
#include <tclap/CmdLine.h>

using namespace GiNaC;
using namespace std;

inline numeric k(numeric s, numeric j_max) {
	numeric res = 0;
	for (numeric i = 1; i <= j_max; i++) {
		res = pow(i, s);
	}
	return res;
}

inline numeric beta(numeric s, numeric j_max) {
	numeric num = 0, dem = 0;
	for (numeric i = 1; i <= j_max; i++) {
		num = num + 1 / pow(i, s-1);
		dem = dem + 1 / pow(i, s);
	}
	return num / dem;
}

inline ex g_z(ex z, numeric s, numeric k_s, numeric j_max) {
	ex res = 0;
	for (unsigned int j = 1; j <= j_max; j++) {
		res+= pow(z, j) / pow(j, s);	
	}
	return res * k(s, j_max);
}

int main(int argc, char **argv) {
	TCLAP::CmdLine cmd("Model Param Simulator", ' ', "0.1");

	TCLAP::ValueArg<unsigned int> n_arg("n", "num", "Set the buffer size", true, 50, "integer >= 0");
  	cmd.add(n_arg);
	
	TCLAP::ValueArg<unsigned int> j_max_arg("j", "j_max", "Set the max numer of packets in a burst", true, 10, "integer > 0");
	cmd.add(j_max_arg);
	
	TCLAP::ValueArg<double> u_arg("u", "u", "Set the fixed rate of the server", true, 10, "double > 0");
	cmd.add(u_arg);
	
	TCLAP::ValueArg<double> lambda_arg("l", "lambda", "Set the average rate of Markov Inter-Arrival times", true, 10, "double > 0");
	cmd.add(lambda_arg);
	
	TCLAP::ValueArg<unsigned int> s_arg("s", "s", "Set the variable s of Zip's law", true, 2, "integer >= 0");
	cmd.add(s_arg);
	
	TCLAP::SwitchArg output_arg("o", "output", "Print output to stdout");
	cmd.add(output_arg);

	cmd.parse(argc, argv);

	unsigned int N = n_arg.getValue();
	numeric j_max = j_max_arg.getValue();
	numeric u = u_arg.getValue();
	numeric lambda = lambda_arg.getValue();	
	numeric	 s = s_arg.getValue();
	numeric ro = beta(s, j_max) * lambda / u;
	
	bool output = output_arg.getValue();

	symbol z("z");

	const ex exp_part = exp( -lambda/u * ( 1 - g_z(z, s, k(s, j_max), j_max)));
	ex p_z = (1 - ro) * ((1 - z) * exp_part) / (exp_part - z);
	
	numeric h = 1E-60;
	for (unsigned int n = 0; n <= N; n++) {
                ex der = 0;
                for (unsigned int i = 0; i <= n; i++) {
                    int sign = i % 2 ? 1 : -1;
                    der += sign * p_z.subs(z == z + (n / 2 - i) * h) / (factorial(n-i) * factorial(i));
		}
		if (output) {
			cout << "[" << n << "]\t" << der.subs(z == 0) << endl << endl;
		}
	}
	
	return 0;
}
